function R = slhitrates(scores, rlabels, qlabels, rankord, rks)
%SLHITRATES Compute the hit-rates in retrieval 
%
% [ Syntax ]
%   - R = slhitrates(scores, rlabels, qlabels, rankord, rks)
%
% [ Arguments ]
%   - scores:       The score matrix (nr x nq)
%   - rlabels:      The labels of references (length nr)
%   - qlabels:      The labels of query samples (length nq)
%   - rankord:      The order of score ranking {'descend' | 'ascend'}
%                   \{:
%                       - descend: rank the references in the descending
%                                  order of score values.
%                                  (larger value indicates better matching)
%                       - ascend:  rank the references in the ascending
%                                  order of score values.
%                                  (smaller value indicates better matching)
%                   \:}
%   - rks:          The ranks at which the hit rates are computed
%                   (rks shoudl be sorted in strictly ascending order)
%   - R:            The resulting hit rates 
%
% [ Description ]
%   - R = slhitrates(scores, rlabels, qlabels, rankord, rks) computes the
%     hit-rates in score-based retrieval. 
%
%     Suppose there are nr references, and nq queries, then the score
%     matrix should be of size nr x nq, in which each column represents 
%     the matching scores of a query with all references. The references
%     can be either samples or a classes/models.
%
%     The retrieval is performed by ranking the references for each query
%     in descending/ascending order of matching scores. 
%
%     For each query, if there exists a reference in the first r references
%     in the resulting ranked list, we say that the retrieval for this
%     query hits at rank r. For a batch of queries, the ratio of the number
%     of queries that hit at rank r to the total number of queries is
%     called the hit rate at rank r. 
%
%     The output R is of the same size as rks, and R(i) equals the hit rate
%     at rank rks(i).
%
% [ Remarks ]
%   - The classification cumulative scores can be computed using this
%     function by inputting rks as a sequence of natural numbers like 1:K.
%
% [ History ]
%   - Created by Dahua Lin, on Jul 19, 2007
%

%% parse and verify input arguments

error(nargchk(4, 5, nargin));

assert(isnumeric(scores) && ndims(scores) == 2, 'sltoolbox:slhitrates:invalidarg', ...
    'The scores should be a 2D numeric matrix.');

assert(isnumeric(rlabels) && isvector(rlabels), 'sltoolbox:slhitrates:invalidarg', ...
    'The rlabels should be a numeric vector.');

assert(isnumeric(qlabels) && isvector(qlabels), 'sltoolbox:slhitrates:invalidarg', ...
    'The qlabels should be a numeric vector.');

[nr, nq] = size(scores);
assert(length(rlabels) == nr && length(qlabels) == nq, 'sltoolbox:slhitrates:sizmismatch', ...
    'The size of scores does not match the number of labels.');

assert(ischar(rankord) && (strcmpi(rankord, 'ascend') || strcmpi(rankord, 'descend')), ...
    'The rankord should be a either ''ascend'' or ''descend''.');

if nargin < 5
    rks = 1 : nr;    
else
    assert(isnumeric(rks) && isvector(rks) && all(rks > 0) && all(rks == fix(rks)), ...
        'sltoolbox:slhitrates:invalidarg', 'rks should be a vector of positive integers.');
    if length(rks) > 1
        assert(all(rks(2:end) > rks(1:end-1)), 'sltoolbox:slhitrates:invalidarg', ...
            'rks should be in strictly increasing order.');
    end
    
    assert(rks(end) <= nr, 'sltoolbox:slhitrates:invalidarg', ...
        'The value of rks should not exceed the number of references');    
end

%% main

% decide size and scheme

rks_start = rks(1);
rks_end = rks(end);
is_rks_sparse = (rks_end - rks_start > 4 * length(rks));

% generate hit-map

[ss, si] = sort(scores, 1, rankord);
clear ss;

if rks_end < nr
    si = si(1:rks_end, :);
end

if ~isequal(size(rlabels), [nr 1])
    rlabels = reshape(rlabels, [nr 1]);
end
if ~isequal(size(qlabels), [1 nq])
    qlabels = reshape(qlabels, [1 nq]);
end

hmap = bsxfun(@eq, rlabels(si), qlabels);
clear si;

% compute hit-rates

if ~is_rks_sparse
    if rks_start == 1
        hitcounts = cumsum(hmap(1:rks_end, :), 1);
    else
        prehits = sum(hmap(1:rks_start-1, :), 1);
        hitcounts = bsxfun(@plus, cumsum(hmap(rks_start:rks_end, :)), prehits);
    end
    
    if length(rks) == rks_start-rks_end+1
        H = hitcounts > 0;
    else
        H = hitcounts(rks - (rks_start-1), :) > 0;
    end
               
else    
    eps = rks(:);
    sps = [1, eps(1:end-1)+1];
    
    hcs = arrayfun(@(sp, ep) sum(hmap(sp:ep, :), 1), sps, eps, 'UniformOutput', false);
    hitcounts = vertcat(hcs{:});
    clear hcs;
    
    H = hitcounts > 0;            
end

R = sum(H, 2) / nq; 
R = reshape(R, size(rks));

